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Abstract 

The Gaussian-filtered Navier-Stokes equations are examined theoretically and a generalized the- 
ory of their numerical stability is proposed. Using the exact expansion series of subfilter-scale 
stresses or integration by parts, the terms describing the interaction between the mean and fluctu- 
ation portions in a statistically steady state are theoretically rewritten into a closed form in terms of 
the known filtered quantities. This process involves high-order derivatives with time-independent 
coefficients. Detailed stability analyses of the closed formulas are presented for determining whether 
a filtered system is numerically stable when finite difference schemes or others are used to solve it. 
It is shown that by the Gaussian filtering operation, second and higher even-order derivatives are 
derived that always exhibit numerical instability in a fixed range of directions; hence, if the filter 
widths are unsuitably large, the filtered Navier-Stokes equations can in certain cases be uncon- 
ditionally unstable even though there is no error in modeling the subfilter-scale stress terms. As 
is proved by a simple example, the essence of the present discussion can be applied to any other 
smooth filters; that is, such a numerical instability problem can arise whenever the dependent 
variables are smoothed out by a filter. 
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I. INTRODUCTION 



Large eddy simulation (LES) 0, HI is one of the typical numerical approaches to tur- 
bulent flows, in which large-scale structures in fluid motion are solved directly while the 
effects of the small-scale eddies are modeled based on a filtering operation that separates 
the high- and low-wavenumber modes in turbulent flow fields. Because LES enables us to 
treat time-dependent, high- Reynolds- number turbulence with substantially smaller compu- 
tational effort and storage than direct numerical simulation (DNS) resolving all scales of 
motion, it has been used in many applications in a variety of research fields, including fluid 
machinery , combustion engineering 0, , atmospheric science 0, 0] , geophysics [H |§, E3| , 
and astrophysics jll. 12l . 

For the past few years, however, several reports have been published that point out the 
incompleteness of current LES modeling. In those repots it has been implied that even 
a completely accurate LES model could be numerically unstable. In Ref. jTj||, Leonard 
showed that the tensor-diffusivity model, re-derived by truncating an exact expansion series 
of subfilter-scale forces [3 EE Q, 

works as a negative-diffusion term in the stretching 
directions of fluid motion and hence, that it could lead to numerical instability when used 
for finite difference schemes. Since that model is exact for first-order velocity fields where the 
velocity components can be described by linear polynomials, the model's negative diffusivity 
can be considered to exhibit the nature of the exact model under a particular condition. 
Also, Winckelmans et al. [TtJ and Kobayashi and Shimomura (3, Tlil li^ pointed out 
that the tensor-diffusivity model behaves unstably in a plane channel flow, also due to the 
model's negative diffusivity. In a comparative study of LES models (the tensor-diffusivity 
and rational LES models) :2lj], Iliescu et al. showed that the tensor-diffusivity model can 
be unstable in a high- Reynolds- number driven cavity flow, as well. In Ref. |22|, Ida and 
Taniguchi derived a closed form of the Gaussian-filtered Navier-Stokes (GFNS) equations 
under a simple assumption about the instantaneous velocity profile and showed theoretically 
that the shears in time-averaged flow fields can also be the seed of the numerical instability 
of the filtered system, because a cross derivative of the filtered velocity component, being 
unconditionally unstable in numerical simulation, appears in the closed formula. The authors 
stressed that the unstable portions in the closed formula must be solved accurately without 
using artificial techniques (e.g., clipping or damping), since those portions derive naturally 
from the filtering operation and are thus a part of the governing equations for LES. In the 
sequel to that paper [23| , Ida and Taniguchi further ascertained that the shears in the time- 
averaged fields can, through the filtering operation, cause the appearance of a numerically 
unstable term that always exhibits a negative diffusivity in a fixed direction, a conclusion 
that is able to explain the problematic instability that has frequently been confronted in 



wall-bounded turbulent flow computations (e.g., Ref. 17, 18]), where a strong shear appears 
in the time-averaged streamwise velocity. The theoretical and numerical findings listed above 
appear to suggest that the filtering operation itself is the underlying cause of the numerical 
instability in LES, raising the question whether a numerically stable LES model can be 
ideally accurate or not. 

A similar scenario can be found in simulation strategies for collisionless plasma kinetics 
that use the Vlasov-Poisson or Vlasov-Maxwell system as a governing equation. In Refs. ^2^, 
25[, Klimas has attempted to apply a Gaussian filter to the Vlasov equations in order to 
mollify the filamentation of the distribution function (an infinitely fine structure in the 
phase space), and found that the filtered Vlasov equations can be rewritten into a closed 
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form in terms of the filtered distribution function and are thus solvable without any empirical 
modeling. (We note here that in Klimas's study the filtering operation was only applied in 
the velocity space, which allows for relatively easy derivation of closed formulas and results 
in only a few additional terms.) In that closed formula, a cross derivative of the filtered 
distribution function appears, which, as Figua et al. suggested |26| (see also Refs. lH 22^). 



makes the filtered system ill-conditioned and unsuitable for numerical simulations using finite 
difference methods or others excluding the spectral method. As with the Navier-Stokes cases 
mentioned above, that finding implies that the Gaussian filtering operation itself, and not 
the modeling or approximations, destabilizes the governing equations. 

The present paper extends the numerical stability analysis of the GFNS equations per- 
formed by Ida and Taniguchi [23| to construct a generalized theory for the numerical insta- 
bility of the system. The present discussion assumes that the flow fields are, as in Ref. |23| . 
in a statistically steady state (an assumption allowing us to decompose the velocity compo- 
nents into time-independent mean and time-dependent fluctuation portions), but the mean 
velocities may be described by high-order polynomials in terms of the spatial coordinates, 
while in Ref. |23| first-order velocity fields have mainly been considered. Under these as- 
sumptions and using the exact expansion series or integration by parts, we rewrite the terms 
that represent the interaction between the mean and fluctuation portions (referred to below 
as "mean-fluctuation terms"), filtered by a Gaussian function, into closed forms involv- 
ing high-order cross derivatives, and show that through Gaussian filtering, various kinds 
of unconditionally unstable terms having time-independent coefficients are derived which 
numerically destabilize the modes in a fixed range of directions. Also, detailed stability 
analyses of the resulting closed formulas are presented to derive a stability criterion for the 
choice of filter widths. In the present paper, for simplicity we only discuss cases in which the 
mean velocity field has one-dimensional (ID) or two-dimensional (2D) structures. Moreover, 
we are not concerned with the commutation error between differentiation and filtering (see 



e.g. Refs. [23, |28|, |29|, |3Cj for recent efforts to resolve the commutation error), assuming 
each filter width to be constant in the corresponding spatial direction. This treatment war- 
rants that the numerically unstable terms that we discuss are not those originating from the 
commutation error, which is a modeling failure. 

The present theoretical investigation has been performed assuming the use of finite dif- 
ference schemes. However, most of the results will be true for other numerical methods (e.g., 
finite volume, finite element, compact differencing) as well. Also, in order to accomplish the 
theoretical investigation without the aid of numerical analysis, the present study neglects 
the cutoff of high-wavenumber modes originating from the use of finite grid spacing. The 
Gaussian filter considered in the present study is, therefore, assumed to approximately rep- 
resent the numerical damping of high-wavenumber modes due to numerical viscosity (also 
originating from the use of finite grid spacing), and also to be an explicit filter applied in- 
dependently of numerical discretization [l7|, |3l|. Because of these assumptions, we use the 
term "subfilter scale" in stead of "subgrid scale" throughout this paper. 

The present paper is organized as follows. In Sec. |HJ the governing equations and defini- 
tions useful for the present investigation are introduced. In Sec. Illl[ the numerical stability 
of arbitrary-order partial differential equations involving high-order cross derivatives is the- 
oretically discussed to derive a stability criterion for them, which is essential for our study. 
Combining the result of this stability analysis with an exact expansion series or integration 
by parts allows us to construct a generalized theory for the filtering instability under sta- 
tistically steady-state conditions. In Sec. II VI several specific examples are investigated to 
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elucidate how the stability criterion restricts the choice of filter widths, and in Sec. El to 
elucidate a fundamental mechanism of the numerically unstable terms and to show that the 
essential part of the present results is true for non-Gaussian smooth filters as well, a simple 
advection problem is considered where the true solution has a discontinuous step. Section 
I VII presents notes on remaining issues that must be resolved to gain a more generalized the- 
ory. As stated in that section, the present theory has several limitations in its applicability. 
We conceive of the present theory as an intermediate step towards a complete theory of 
the numerical instability of the GFNS equations. Section I VI II summarizes this paper, and 
the appendix presents a mathematical proof of the exact expansion series using elementary 
mathematics, thereby assuring the self-consistency of the present paper. 



II. GOVERNING EQUATIONS, FILTERING OPERATIONS, AND DEFINI- 
TIONS 

Incompressible viscous fluid flows are described by the Navier-Stokes equations: 

duj dujUj = _dp_ ^ d 2 Uj for ^ = x 2 3 (1) 
dt dxj dxi dxjdxj 

where the summation convention is assumed, and Ui (i — 1, 2, 3) are the velocity components, 
p is the pressure divided by the constant fluid density, and v is the kinematic viscosity. The 
subfilter-scale terms, resulting from a low-pass filtering operation, are derived from the 
convection terms (i.e., the second term of Eq. (|T|)). In what follows, we assume that the 
velocity components can be decomposed into time-averaged and fluctuation portions as 

u(x,t) = U(x) +u'(x,t), (3) 

where u = (ui, ^2,^3), U = (^1,^2,^3), and u' = (u[, u' 2 , u' 3 ) . We also assume that the 
mean velocity U is time-independent; i.e., that the flow is in a statistically steady state. 
From Eqs. (J2J) and 0, one can derive 

OXi OXi 

Using Eqs. (j3J) and (jlj), the convection term in Eq. (JTJ) is rewritten as follows: 

dujUi , du- dUi 

— = h i + u j — - + Uj — , 5 

OXj J OXj OXj 

h^U^ + ^u'j fori = 1,2, 3, (6) 

where hi represents the interaction between the mean and fluctuation portions, on which we 
focus our attention. 

The filter function is assumed to be Gaussian: 



7rA 2 ^ V A 2 
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which satisfies J^~^° Q(X;A)dX = 1, where 7 is commonly set to 6 in LES and A is 
the filter width. Using this, the filtering operation in the X{ direction is performed as a 
convolution integral: 



X=oo 

f{x h . . . , t) = J g( Xl - X; Ai)f(X, t)dX 

X=~oo 

= Gi*f, 

where the overbar denotes the filtered quantities. Three-dimensional (3D) filtering is 
achieved by successively performing this convolution as follows: 

X=oo 3 

/(x, *) = / n Gfa - Ai)/(x, t)dx = & * [g 2 * (g 3 * /)] 

X=-oo i=l 

= Gl23 * /, 

where Aj (i = 1, 2, 3) is the filter width in the Xi direction. As stated in Sec. |U we assume 
throughout this paper that each filter width is a constant, and thus 



df d(Gi*f) 



Gi * {Gj * /) = Gj * {Gi * /) for i,j = 1, 2, 3. 

For the convenience of the following discussion, we introduce the residual stress function 
defined as 



n a [F(y 1 ,u' 1 ,...)] = g a *F(u 1 ,u' 1 , 

which yields, for example, 



-F(0 o *£/i,0 o *ui,...) 



(7) 



Hi 



3 8Xn 



Gi*[U. 



dxj 



- (Gi * Uj 



d(Gi*u2 
dx; 



(8) 



Tll23 [UW2} = Gl2-i * {.UW2) ~ (^123 * «'l)(£l23 * «a) 



(9) 



Based on the above assumptions and definitions, the Navier-Stokes equations filtered 
using a 3D Gaussian filter are written as 



du[ diijUi 
dt 8xa 



dp d 2 u 
+ v 



dx 



n 



123 



dujUi 
dxj j 



(10) 



or 



du[ dujUi 

~dt + ~~~~~ 



dXj 



dp <9 2 w' 
+ v 



dx. 



dxjdxj 



Tl\2z [hi 



K 



123 



3 dxj 



+ v 



d 2 U t 
dxjdxj 



71 



123 



3 dxi 



(11) 



which can be considered an equation for both and Ui because du^/dt = dui/dt. The 
terms in the last parentheses of Eq. (fTTj) can be considered time-independent source terms, 
which may have no influence on the numerical stability. The next to last term represents 
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the residual stress forces due to the nonlinear interaction between the fluctuation portions, 
the stability analysis of which is difficult to complete theoretically and thus requires nu- 
merical experiments. Although it has been pointed out that terms having the same form as 
'R<V2s[ v >' : jdu' i / dxj] can instantaneously be numerically unstable [HI sucn terms should not 
necessarily lead to numerical instability in actual computation, because their time-averaged 
nature can be dissipative. In what follows, we only consider the numerical stability of 

d 2 u' 

"s^ (12) 

i.e., the difference between the molecular viscosity and the residual stress forces due to the 
mean-fluctuation interaction, and do not take into consideration the numerical effects of 
the nonlinear term. In this respect, our theoretical investigation is incomplete. Comments 
on potential approaches to resolving this incompleteness are given in Sec. As shown in 
what follows, the closed form of Eq. (|12|) has time-independent coefficients, meaning that 
the numerical stability of this portion is time-independent. 

We introduce here some mathematical tools that allow us to rewrite the filtered mean- 
fluctuation term 7^123 [hi] into a closed form. Yeo [li| and others [l3|, EH 3] have derived 
a very interesting identity, which is applicable to all differentiable and continuous functions 
f(x) and g(x), 

(f9)-f9 = J2^\-^) d^Q^- ( 13 ) 

n=l v ' 

Here, the overbar indicates rr-directional Gaussian filtering. It is worth noting that the 
right-hand side of this identity only involves the known filtered quantitie s / a nd g. This 
outstanding feature of the series allows us to rewrite the residual TZ[fg] = (fg) — fg into a 
closed form. For 2D Gaussian filtering in the (xi,X2) plane, this series becomes 

t=i V 2 7 J Sx k dx k ^ 2! V 2 7 J \ 27 J dx k dx, dx k dx, 

+ YL(£)(*i)(*l)-*l *!— + .... (i4, 

k f^ =1 3- V 2 7 / \27/ V 27 / dx k dxidx m dx k dxidx m 
The following identity is also useful: 

Qi * (xj) = (xi + ^-^j (ft * /), i = 1, 2, 3, (15) 



which can be derived using integration by parts (e.g., Refs. ji^, 2^, 2$, 32|). This yields, for 
example, 

Gi-kXi = Xi, (16) 

ft*^ = ^ + ^, (17) 
27 

3A 2 

Qi*x\ = x\ + — -Xi, (18) 
27 



6 



where the summation convention is not adopted. Using Eq. (|TE|). Eq. (|T5|l can be rewritten 
into 

*^ = ^r- i = 1 ' 2 ' 3 ' (19) 

The expansion series (|T3j) and (JT4*j) enable us to derive a closed form of Eq. (JT2*j) . As can 
be seen from these series, the closed form has high-order derivatives including high-order 
cross derivatives, and the coefficients of these derivatives are time-independent, such as Uj 
and dUi/dxj in hi. The numerical stability of such high-order derivatives are examined 
below. 



III. NUMERICAL STABILITY OF ARBITRARY-ORDER PARTIAL DIFFEREN- 
TIAL EQUATIONS 

We derive and examine in this section a numerical-stability criterion to determine whether 
an arbitrary-order partial differential equation (PDE) can be solved stably (and accurately) 
by a finite difference scheme. Although the numerical stability of PDEs is known to depend 
on the discretization scheme applied, we do not discuss a certain form of finite differencing. 
We instead consider the exact amplification factor of the PDEs, which is essential and may be 
sufficient for our aim. It is well known that a diffusion equation, for example, is numerically 
stable (numerically unstable) when the coefficient of the diffusion term is positive (negative), 
i.e., when the exact amplification factor is less than (greater than) unity. We assume here 
that such is also the case for other types of PDEs including high-order ones, and use their 
amplification factors to judge whether a stable finite difference scheme can exist for the 
corresponding PDE. 

Let us consider the exact solution of a 2D arbitrary-order PDE, 

s-'(k)"K)"'- <w 

/ = f(x,y,t) 

where fi is a real constant and m, n > are integers. The initial condition is 

f(x,y,0) = exp(i(k x x + k y y)), (21) 

where k x and k y are real constants, i is the imaginary unit, and the boundary conditions are 
periodic. Suppose that the exact solution of this PDE has the form of 

f(x,y,t) = eMut)f(x,y,0), (22) 

where u is a complex constant. Substituting this into Eq. (|2U|) yields 

u = i™+>£;™A£. (23) 

The characteristic of this exact solution can roughly be categorized into the following two 
solutions: 

For m + n = 2M: 

f = exp((-l) JW 'nh™h»t) exp{i{k x x + k y y)), (24) 
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For m + n = 2M + 1: 



/ = exp(i(-l) M /i^^) exp(i(^x + kyy)), (25) 

where M = 0, 1, 2, • • • . The former represents exponential decay or growth of the solution, 
whereas the latter represents phase shift without changing amplitude. 

Equation (|24j) suggests that in order for a stable finite difference scheme for m + n = 2M 
to exist, the amplification exponent of solution (I24J1 . 

a(m } n^) = (-l) M ^k n y , (26) 

must be zero or less for any value of wavenumbers. Furthermore, if m and n are even 
numbers, k™k y in this exponent is always positive and the stability is thus determined by 
the sign of (— l) A/ /i. (If, for example, m = 2 and n = 0, which results in M = 1, the PDEs 
with \i > are numerically stable, while those with fi < are unconditionally unstable, 
a conclusion that is consistent with the well-known fact that a positive diffusion equation 
can be solved stably but a negative one can not.) Otherwise, i.e., if m and n are odd 
numbers, k™ky can be either positive or negative, meaning that the PDEs in this case are 
always unconditionally unstable because modes in any direction can appear in turbulent 
flows. (This conclusion is consistent with the known fact that the PDE for m = n = 1 is 
unconditionally unstable; see, e.g., Refs. [HI Hfjf). For (— fi < 0, for example, the 
modes of k x k y < must be unstable, implying that if the sign of \x is constant, the PDEs 
in this case always exhibit numerical instability in a fixed range of directions. This finding 
plays an important role in our main subject discussed in the next section. 

On the other hand, Eq. (|23j) suggests that if m + n is an odd number, a stable finite 
difference scheme should always exist because the amplification exponent 

a(m,n;/i) = \{-\) M ^k n y (27) 

has an imaginary value and hence the absolute value of the amplification factor, 
| exp(i(— 1) M fik™kyt)\, is unity. Indeed, the advection equations (corresponding to the case 
of, e.g., m = 1 with n = 0) and the Korteweg-de Vries (KdV) equations involving third- 
and/or fifth-order dispersion terms (e.g., for m — 1, 3, 5 with n = 0) have been solved stably 
and accurately by finite difference schemes; see, e.g., Refs. @, H @ for recent progress m 
finite difference schemes for KdV equations. 

The present theoretical results can be summarized as follows: A stable finite difference 
solver must exist for odd-order PDEs (i.e., when m + n is an odd number). For even-order 
PDEs, a stable solver exists only if both m and n are even numbers and (— l)( m + n )/ 2 ^ is neg- 
ative; otherwise, the PDE is unconditionally unstable by any finite difference scheme, since 
the numerical perturbations grow exponentially in numerical simulations. This conclusion 
may also be true for the finite volume, finite element, and compact difference schemes. 

The total amplification exponent «t of a complicated PDE, 

can be determined by 

«t = ^ Re[a(m i; ^)]. (29) 

i 

Though this is, for variable //j, only an approximation, it should work sufficiently in many 
situations. 



IV. STABILITY ANALYSIS OF THE FILTERED SYSTEM 



We present several analytical results for the numerical stability of a filtered system deter- 
mined by combining the stability analysis described in the previous section with the exact 
expansion series (|13|) and (|14j) or the identity given using integration by parts, (|15j) . As 
stated in Sec. |H1 we consider up to 2D cases for the sake of simplicity, and only analyze the 
stability of Eq. (fT21). 



A. ID mean velocity cases 

Suppose that 

£/i = £/i(x 2 ) and U 2 = U 3 = 0, 
(this means that the velocity satisfies the divergence-free condition), which leads to 

< 3 °> 

h 2 = (3D 

*. = U*. (32) 

Because U\ depends only on x 2 , filtering in the x\ and x 3 directions (i.e., in the homogeneous 
directions) results in 

Ki 3 [hi) = fori = 1,2, 3, (33) 

When the Gaussian filter in the x 2 direction is applied to Eqs. (|3(Jj) - (j32j) . some mathemat- 
ical manipulations are needed to obtain a closed formula because U\ and dUi/dx2 cannot 
simply be put outside the convolution operation. We consider here the case where U\ is 
described by a finite-order polynomial: 

TV 
n=0 

where a n (n — 0,1, ... , N) are real constants and N is the order of this polynomial. Using 
the ID expansion series (|13|). we have 



where 



TZ 2 [h 2 ] = C[i]u' 2 , (36) 

Knifa] = >C[i]«3, (37) 
/ A|\ n d n Ui 

[1] ^nl^J dx n 2 dx x dxf 1 ' 
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°° 1 /A 2 \ n d n+1 Ui d n 
C[2] ^^[n\\^) ~d^ r d^' (39) 

In the following, we examine some special cases for deriving the stability conditions for the 
ID flows. 

If N = 1, operators (JHHj) and (|3*9]). respectively, reduce to 

c AA|\gg, d 2 
^' \27/ dx 2 dx\dx 2 



C[2) = 0, (41) 



and Eq. (|T2j) then becomes 

d 2 



c ■)"'- 1 L±:l (42) 

Here the diffusion operator in the X3 direction is neglected because it does not alter the 
resulting stability condition that is applicable to any wavenumber; note that the neglected 
operator does not stabilize the modes in the (07,22) plane but that C[i], being unstable, 
only has derivatives with respect to x\ and x 2 . Based on Eqs. (j2fij) and (|29j) . we have 

a T ~ -u(kl + kl) + (^P) ^k x k y . (43) 



Substituting (k x ,k y ) = |k| (cos 9, sin 8) with G [— ir, it] into ay < yields 



A'ju dx 2 J 471/ 



2 



, A 2 dUt A A 2 
1 > max —2-— A sin 20 ' 



(44) 



This stability condition is equivalent to that for linear shears determined in Ref. [23j by a 
different approach, which imposed a strong restriction on the choice of the wall-normal filter 
width for use in the viscous sublayer in plane channel flows. 
For N = 2, £[1] and C[ 2 ] become 

^A 2 \ dUi d 2 l/A 2 \ 2 d 2 t/i <9 3 



^ \27/ dx 2 dxidx 2 2 \ 27 / cfeicte 2 .' 
/A 2 \ 9 2 t/i d 



^=\^)m^- (46) 

As proven in Sec. 11111 the third-order operator in Eq. (J45|) can be ignored in stability analysis. 
Moreover, C[ 2 \u' 2 in Eq. can also be neglected, because the numerical stability of u 2 is 
determined independently of Eq. (JHSJ), by Eq. (jSEJ), and furthermore, if Eq. (jHEJ) is unstable, 
then Eq. (|35p should be unstable. Therefore, the stability condition in the present example 
is the same as that in the previous case, Eq. (jH)) . 

For N = 3, the total amplification exponent for v[d 2 jdx\ + d 2 jdx\) — £m is 

- - -*« + « + (g ) gw, - i (f ) 3 gwj- (47) 
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Using this and (k x , k y ) = |k| (cos9,sm9), we obtain the following stability condition: 

1 [A 2 \dU x 1 /A 2 \ 3 d^Ux „ l2 , , — . , . 

1 > — — — -s — — 4- k s i ± Vl - s 2 , 48 

-2u\2-fjdx 2 24z/ V 27 / dx\ 1 1 V ; ' V ; 

s = sin 2^ G [-1,1]. 

This has to be fulfilled for any choice of |k| and 9. 

To show how the restriction (|4~H|) works in a realistic situation, we consider here the 
inertial sublayer forming in a plane channel flow. Following Dean |36j , the streamwise mean 
velocity in the inertial sublayer is approximately described by 

= 2.44 ln(x+) + 5.17, (49) 

where u T is the wall-friction velocity and x 2 = X2U T / v is the distance from the plane wall in 
wall units. For 30 < x 2 < 80, Eq. (|4~9*|) can be well approximated by a cubic polynomial: 

^Ml = j-a n xt\ (50) 

T n=0 

a = 10.5, ai = 0.134, 

a 2 = -1.23 x 10" 3 , a 3 = 5 x 10~ 6 . 

Using this and Eqs. (fT6|) and (fTYj) . the filtered derivatives in Eq. (}48|) are determined as 

(51) 



dUi u T 

a — = — 

OXo V 



A +2 

ai + 2a 2 x 2 + 3a 3 I x 2 2 + 



d 3 ^ 



Ut 



2 1 



6a 3 , (52) 



(53) 



where A^ = A 2 u T / v. Substituting them and 7 = 6 into Eq. (jjHj) yields 

/AM1A+ 2 a 3 A+ 4 \ a 3 A+ 6 ., l2 /u r \-2 , , 

1 > [ 1 2 + 6 2 s - 2 k ( — ) s 1± Vl-s 2 , 
- V 24 96 / 6912 1 U/ v ; ' 

where 

A[i] = ai + 202^^ + 3^0;^ 2 . 

Assuming that max(|k|)A 2 = vr, i.e., the maximum resolved wavenumber is determined by 
the Nyquist wavenumber based on the wall-normal filter width, Eq. (JSlIj) can be further 
rewritten as 

/A [11 A+ 2 a 3 A+ 4 \ a 3 7r 2 A+ 4 , , , , 

1 > + s - — — s l±Vl-s 2 . 54 

- V 24 96 J 6912 v ; v ; 

At x 2 = 55, for instance, this becomes 



1 > (1.85 x 10~ 3 A +2 + 5.21 x 10~ 8 A +4 )s - 7.14 x 10~ 9 A +4 s(l ± Vl - s 2 ). (55) 
From this, for A + = 10, 20, and 25, respectively, we have 



1 > 0.185s - 7.14 x 10~ 5 s(l ± Vl - s 2 ), (56) 
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1 > 0.748s - 0.00114s(l ± Vl-s 2 ), 



(57) 



1 > 1.18s - 0.00279s(l ± Vl - s 2 ). (58) 

The first two are true for any s G [—1,1], but the last is not. (Note that max s [s(l ± 
\J\ — s 2 )] = — min s [s(l ± \/l — s 2 )] = 3\/3/4 ~ 1.30.) The filter width suggested here for 
stability is comparable to that used in actual channel flow computations. 

Based on the stability analysis described in Sec. IIIII it is found that for larger N, all of the 
even-order differential operators in Eq. (|3*Hj) are unstable, whereas the odd-order ones have 
no influence on the numerical stability. That is, the high-order terms do not help stability. 
This result suggests that in most cases of ID shear, the subfilter-scale stress terms would 
be unstable, thus leading to a divergence of numerical solution, if an unsuitably large filter 
width is used. 



B. 2D mean velocity cases 

Next, we consider 2D problems. Suppose that 

U 1 = Ui(x 1 ,x 2 ), U 2 = U 2 (x l ,x 2 ), and U 3 

resulting in 



h = Vu[ + — l -u\ + -M, (59) 



dx i dx 



2 



and 



dxi dx\ dx 
Here we introduced an advection operator, 

d d 
V = U 1 — + U 2 - 



h = Vu' 3 , (61) 

dUi dUi dU 2 n 

+ — 2 - = 0. (62) 



dxi dx 2 

Because the mean velocity is independent of x 3 , we know that 

K 3 [hi) = fori = 1,2, 3, (63) 

and consequently the resulting formulas of the residual stresses for £123* and for Q\ 2 -k are 
the same, allowing for the consideration based on 2D filtering. However, even in 2D, the 
complete set of the closed residual forces derived using Eq. ()14J) is intricate and inconvenient 
for theoretical analysis, and hence we only consider some simple cases. 

The first example assumes that the mean velocities are described locally by 

U\ = bxi and U 2 = —bx 2 , (64) 
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where b is a positive constant, that is, stretches uniformly in the x\ direction. Here, the term 
"locally" means "in a region sufficiently larger than the filter widths." These assumptions 
reduce Eqs. ff5T3J)-(f6"Tl) to 

h = Vu[ + bu[, (65) 



h 2 = Vu 2 — bu' 2 , 



(66) 



Since, in this case, 



Gi * U, 



dx k 



h-* = T>u'n. 



U k a{g '*<> for^k, no summation over * 

OXk 



(67) 



and 



■R 12 [bu\] =K 12 [bu' 2 ] =0 
are true, the residual forces are expressed as 



bx\ 



bx- 



dxo 



1,2,3. 



(68) 



Using the ID expansion series (J13)) or integration by parts (|15|) . Eq. (|68|) is rewritten into 
the closed form, 

Aj d 2 u'i A 2 2 d 2 ^ 



Tti 2 \hi 



-b 



27 dx\ 



(69) 



27 dx\ *-j ^ 2 

This acts as negative diffusion in the x\ direction (i.e., in the stretching direction) but 
as positive diffusion in the x 2 direction, a result that is consistent with Leonard's finding 
shown in the studies on tensor-diffusivity models 0, 0|- The amplification exponent of 
the difference between the viscosity term and 1Zi 2 [hj\ is 



v b^\ 



v + b^ 2 ) ir 

27, 



which expression leads to the stability condition 

x/27^ 



Ai< 



Vb • 



(70) 



(71) 



Equation (fTTj) indicates that a smaller filter width is needed for stronger stretching, and the 
largest filter width usable in the stretching direction is inversely proportional to the square 
root of the velocity gradient. 

The next example is complicated; not only the normal stresses but also a shear stress 
appear in the mean field. Suppose that U\ is described locally by 



U 1 (xx,x 2 ) = bxix 2 , 



(72) 



which involves both normal and shear stresses. Then one has, from the divergence-free 
condition, 

b o 



U 2 (x 2 ) 



-.x 



2- 



(73) 
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where U 2 (0) = is assumed without loss of generality, and b denotes a positive constant 
as in the previous example; see Fig. ^ showing the vector plot of this velocity field around 
the origin (xi,x 2 ) = (0,0). The upper side (x 2 > 0) of this figure seems to represent a flow 
impinging on the wall located at x 2 = 0. Substituting Eqs. (f72|) and (f75|) . Eqs. (J5HI)-(JHD) 
become 

hi = Vu\ + bx 2 u[ + bxiu' 2 , (74) 



h 2 = Vu 2 — bx 2 u 2 , 



(75) 



Vu' 3 . 



(76) 



Since 



G12 * Ut 



ok 

dxi 



bGi* 



X1G2 * x 2 



Eq. (|T5)l (and also can be used to obtain 



K 



12 



dxi 



-bx 



+ 



'-bxx 



' dxi 



d 2 u> , A 2 , A? «9 3 ^ 



+ 



27 dx\ 27 dx\dx 2 27 27 dx\dx 2 



(77) 



where we used £/i 2 •kU\ = U\. Furthermore, using Eq. (fE?j) yields 



ft 



12 



?7< 



dx 2 



U, 



dx 2 



A 2 d 2 v! t 1 /A 2 ^ 2 ^ 

-bx 2 ^-^r — b—^. 

2 7 dx\ 2 V 27 J dxi 



where we used 



dU 2 dU 2 



dx 2 dxo 



-bx? and 



d 2 U 2 
dxi 



d 2 u 2 

dxi 



-b. 



(78) 



(79) 



-2 u ^2 

The remaining terms can also be rewritten into a closed form using Eq. (j!5|) or (j!3|) . Finally, 
we obtain the closed residual stresses: 



^ 2 ^^^ + ^^ + 6 2^^ 

. An dw'o 



7^X2 [h 2 ] = £[3}U 2 

n 12 [h 3 _ 



27 dx 2 

3] "3) 



(80) 

(81) 
(82) 



where 



[3] 



6X5 



+ 



<9 2 



A|^ 2 __ A|^ 

27 <9x 2 27 dx\) ' 27" J dx\dx 2 



A 2 



A 2 /A 2 a 3 1A 2 <9 3 
H -& ' 1 



27 \ 27 dx\dx 2 2 27 dxi 



(83) 



In £[ 3 ] we can see various kinds of operators: negative and positive diffusions, second- and 
third-order cross derivatives, and third-order dispersion, among which third-order terms do 
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not alter the amplification exponent olt- Also, the last two terms of Eq. (j80j) and the last 
of Eq. (|81jh first-order derivatives, may not concern the numerical stability. That is, the 
differential operators responsible for stability are 



v 



d 2 



+ 



d 2 



dx\ dx 



— bx- 



A 2 d 2 

27 dx\ 27 dx\) 27 ™ y dx\dx 2 



A 2 d 2 



A 2 d 2 
-bx\- 



whose amplification exponent is 



From this, the stability condition of the present example is determined as 



(84) 



-v{k 2 x + k 2 y ) + bx 2 (J^k 2 x - ^k^j + ^bx x k x k y . (85) 



1 > -^—[{A\ - A^)x 2 + (Al + Al)x 2 cos 29 + A\xi sin 29], 



(86) 



which should be fulfilled for any choice of 9. 

Let us consider some particular cases to show how Eq. (|86|) works. The 2D formula can 
be used to discover the stability conditions for ID filtering as well, because lim {Qi* f) = f. 

For A 2 — > and A x — > 0, respectively, Eq. (jHo^) reduces to 



1 > - — {x 2 + x 2 cos29)Al, 
47 z/ 



(87) 



1 > 



47Z/ 



'— x 2 + x 2 cos 29 + x\ sin 29) A 2 ,. 



(88) 



On the other hand, if the condition Ai = A 2 = A needs to be satisfied for some factor, then 
Eq. becomes 

cos 29 + Xi sin 29) A . 

A^v 

Below we briefly discuss these three cases. 
For Xl = x 2 > 0, Eqs. (HZI)-© reduce to 



1 > max 



1 > max 

bx 2 ^ 
A^v 



bx-) 



(1 + cos 29) A\ 



A^v 

1 + cos 29 + sin 29) A 



and 



1 > max 



bx 2 
Ajv 



[2 cos 29 + sin 29) A 2 



bx 2 
A^v 

bx 2 
A^v 

bx 2 
A^v 



2A?, 
\ (v^-l)A 
v^A 2 , 



and the respective influential differential operators are 

a\ o 2 



v 



d?_ d_i 

dx\ dx\ 



bxn 



27 dx\ 



2 ■ 



(89) 



(90) 
(91) 
(92) 

(93) 
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(&_ _cf_\ bx _ a 2 a 2 

y&r 2 cfe 2 , / 2 27 ftr 2 . 2 27 dxidx 2 ' 

^ i/f— — ^ -6a; f— — -— — ^ -bx — — — (95) 

\dxf dx\) 2 \2^dx\ 27 dx\) 2 27 dx\dx 2 ' ' 

Among the stability conditions, Eq. (J9l"| gives the weakest restriction on the filter width; the 
second-to-last term of Eq. (|94)l . being a positive diffusion term resulting from compression 
in the x 2 direction, mitigates the instability of the last term, the cross derivative resulting 
from the shear in the same direction. In contrast, Eq. (J92)) is the most restrictive condition, 
resulting from the coexistence of a negative-diffusion term and a cross derivative term. 
For x\ = x 2 < 0, on the other hand, Eqs. (|87j) - (|89j) become 

1 > x A 2 , (96) 



+ (97) 

1 >- ^ w 

In this case, Eqs. (|97|) and ()98|). whose respective influential derivatives have both negative- 
diffusion and cross-derivative terms, give restrictions of almost equal strength, while Eq. (|96jl. 
involving positive-diffusion terms only, imposes no restriction. The results provided here 
denote that the numerical stability of the filtered system and the possible choice of filter 
widths depend on how the filtering operations are applied; this conclusion confirms the same 
assertion presented in Ref. [23]. 



V. DISCUSSION OF THE NEGATIVE-DIFFUSION TERM 

As has been shown in the previous sections, many kinds of numerically unstable terms 
are derived by the Gaussian filtering operation. In this section we would like to remark 
on the negative-diffusion term appearing in the stretching direction to clarify why such 
unstable terms appear and how the terms work. The discussion also clarifies that the 
present su gges tions are basically true even for other smooth filters. 

In Ref. |l3j, Leonard considered the pure advection of a sinusoidal wave in a stretching 
velocity field where the amplitude of the (unfiltered) sinusoidal wave remains constant, 
and provided an interpretation of the negative diffusivity of the tensor-diffusivity model as 
follows: As a sinusoidal wave propagates into a stretching velocity field, its wavenumber 
gradually decreases, resulting in the increase of the Gaussian-filtered value of the amplitude 
because of the larger value of the Gaussian filter function for a lower wavenumber. The 
negative-diffusion term represents this amplification of the filtered value resulting from the 
wavenumber shift. We introduce here a different interpretation of the negative diffusivity, 
using Fig. El Let us consider the ID pure advection of a step function, 



f(x,t = 0) 



1 for x < Xq, 
otherwise, 
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in a stretching velocity field u = ax, where xq (> 0) is the initial position of the discontinuity 
in the step function and a is a positive constant. The exact solution of df /dt + udf jdx = 
under this condition is 



which indicates that the discontinuous step will change its position without changing its 
height and profile; that is, the wavenumber shift does not occur in the unfiltered true 
solution of the present example. Applying the Gaussian filter to this solution smoothes out 
the discontinuity to yield a mollified step whose characteristic width is about 2 A, where A 
is the characteristic length of the applied Gaussian filter. Obviously the characteristic width 
of the mollified step is time- independent if A is constant. If, however, the pure advection 
equation in terms of the filtered value, df/dt + udf/dx = 0, is used to advance the filtered 
profile (corresponding to the case where the residual stress term is clipped) , the width of the 
mollified step, unlike that in the true solution, increases gradually as time goes by due to the 
stretching velocity field where a downstream fluid particle moves faster than an upstream 
particle. To counteract this artificial expansion of the mollified step, a modification by 
negative, not positive, diffusion is necessary, which sharpens /. In the case of a compression 
velocity field, a similar but opposite treatment, i.e., the addition of a positive diffusion 
term, is needed because an artificial compression of the mollified step arises if only the pure 
advection equation is assumed. 

The present physical picture may allow us to conclude that the subfilter-scale terms 
should have an analogous negative diffusivity also for any other filter functions that smooth 
the profile of dependent variables. In the above discussion, there is no reason that the filter 
shape must be Gaussian. The artificial expansion of the step discussed above must occur 
whenever the step is smoothed out by a smooth filter but the pure advection equation is 
solved. The above discussion also suggests that a numerically unstable term can appear by 
filtering even if the true solution (both filtered and unfiltered) is physically bounded. 

VI. TOWARDS A FURTHER GENERALIZATION 

The present theory has several limitations in its applicability resulting from the assump- 
tions and simplifications made. In this section we remark on some significant issues that 
must be resolved for further generalization. 

To construct a more general theory for the numerical stability of the GFNS equations in 
statistically steady states, one has to elucidate the numerical stability of, not only Eq. (|2*Sjl . 
but also 



coefficients, determined by the expansion series, and in general 

C[ ab ]C[ cd ] ^ C[ cd ]C[ab] for (a, b) ^ (c,d); 

that is, these operators do not commute with each other. Equation (|99|) represents a com- 
plicated coupling between the equations of different velocity components. In the present 
study, having assumed ID or 2D mean velocity fields, some of the mean-fluctuation terms 
were uncoupled as Eqs. (0TJ), (|52*|) . (pT|). and (JZBJ), and hence the determined stability con- 
ditions are accurate only for the corresponding uncoupled portions. In fully 3D cases where 
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high-order velocity fields must be assumed, however, we may well have to treat the fully 
coupled system (J^H)). 

The strongest assumption among those made in this paper may be the omission of the 
nonlinear fluctuation term 7^123 [u'jdu'Jdxj] in Eq. (JTTJ). We could not take into consideration 
the effects of this term since the theoretical investigation of it is quite difficult to perform 
accurately, and the present theory is thus only valid when the fluctuation is small. Terms 
of this type, as is well known, have a dissipative character in many turbulent flows, and 
hence when the dissipation of the omitted term is strong enough, the instability of the 
mean-fluctuation terms can be eliminated completely. One possible way to gain detailed 
knowledge of the nonlinear term is an a prioritest using DNS, which enables us to determine 
the values of all terms in the GFNS equations. Observing and examining the numerically 
determined terms should allow us to obtain a more accurate prediction of the numerical 
instability. When, for example, the absolute value of the nonlinear term (plus the molecular 
viscosity) is smaller than that of the sum of the unstable terms, the instability can not be 
eliminated irrespective of the specific characteristic of the nonlinear term. Also, comparing 
the amounts of the energy dissipations due to the unstable terms and the nonlinear terms 
should provide a useful insight. This issue will be addressed in a future paper. 

Lastly, we make a brief comment on cases with a nonuniform filter width. Consider 
again pure advection of a step function being smoothed by a smooth filter. When the filter 
width is spatially nonuniform, in the advection process the width of the mollified step varies 
according to its position even for a constant velocity, and hence a negative diffusion must 
take place at least in the period where the width decreases. Such an effect of nonuniform 
filtering must be discussed carefully in the near future. 

VII. CONCLUSION 

We have presented a generalized theory for the numerical instability of the Gaussian- 
filtered Navier-Stokes equations. The theory allows for high-order mean velocity fields and 
high-order derivatives resulting from the Gaussian filtering operation. Also, we have de- 
scribed stability conditions regarding the choice of the filter widths in several situations, the 
violation of which should lead to unconditional numerical instability of the filtered system 
even when a completely accurate subfilter-scale model exists and is used. It is worth noting 
again that the closed formulas of the filtered mean-fluctuation terms determined under sta- 
tistically steady-state conditions involve various kinds of unstable derivatives that, because 
their coefficients are time- independent, always exhibit numerical instability in a fixed range 
of directions. As has been proven by a simple example, the essential part of the present 
results can be true even if a non-Gaussian smooth filter is assumed. 

We stress that if one skirts this numerical instability problem, the accuracy of the LES 
results will plateau. It is hard to imagine that ideally accurate solutions can be achieved by 
incorporating an artificial damping or clipping technique to avoid this numerical difficulty, 
because when the subfilter-scale terms act unstably, the absolute values of their unstable 
portions must be greater than that of the molecular and turbulent viscosities, and the 
adoption of such artificial techniques thus corresponds to the disregard of a term whose 
dominance is greater than that of a term involved ab initio in the Navier-Stokes equations. 
Recently, Moeleker and Leonard 0] have tackled this numerical instability problem and 
proposed an approach to potentially resolve it, based on an anisotropic particle method 
incorporating a remeshing technique. Their method has provided excellent results for a 
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2D scalar advection-diffusion equation with a known velocity field. However, the extension 
of that approach to the 3D Navier-Stokes equations has to the author's knowledge not 
yet been achieved. Because finite difference schemes have been used widely in turbulence 
computations, constructing a stable and accurate solver in the finite difference framework 
would be preferable, though it will be an exceedingly difficult task and might even be an 
unsolvable problem, such as the gravitational three-body problem and the algebraic solution 
of general fifth-order polynomial equations of one variable. We do not know so far whether 
this instability problem is resolvable or not, but we can say that this problem is not something 
that can be avoided when an accurate solution is desired. 
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APPENDIX: ALTERNATIVE DERIVATION OF THE EXACT EXPANSION SE- 
RIES FOR GAUSSIAN FILTERS 

The exact expansion series for Gaussian filters has served as a powerful tool in our study. 
We present here a derivation of the series to improve the self-consistency of the present paper. 
This derivation seems to be rather intricate and drawn out compared to those by Moeleker 
and Leonard 0] and by Carati et al. [15], but it only consists of elementary mathematics: 
the Taylor expansion, integration by parts, and some simple algebraic operations. Some 
readers may prefer the present derivation. 

Let a(x), b(x), f(x), and g(x) be arbitrary, different iable and continuous functions of x. 
For the Gaussian filter with 7 = 1/2 and the characteristic width of A, the exact expansion 
series in ID reads 



In what follows, we derive the right-hand side of this equation from the left-hand side. 
Taylor expanding a with respect to x results in 




(A.l) 



n=0 




oc 



ml dx' 



(A.2) 



m=0 



2=0 



Substituting this into (ab) yields 




(A.3) 



Successively using 



(*/) 



Cf with C = x + A 2 — , 



(A.4) 
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which corresponds to Eq. (|T5j) derived using integration by parts, we obtain the identity [2^ 



(x™f) = C m f. (A.5) 

This rewrites Eq. (jA.3|) as 

oo 

p) = ^[A m (£ m 6)]. (A.6) 



m=0 



The component £ m 6 in Eq. (jA.6|) can further be rewritten as follows: Operating C once 
on h yields 

Cb = xb + A 2 -^-b 

ox 

. (£ M,(A^)V(£M)(A^)\ (A.7, 

where 

C l • 1 = (x + A 2 -^] l = x and £° • 1 = 1. (A.8) 



<9x 

We introduce here 

Bn,m = (£ iV ■ 1) ( A 2_ ht ) 6. (A.9) 
Based on Eq. (jA.7j) . definition (jA.flj) . and 



C{sf) = (Cg)f + g[A^)f 1 



the following identities are derived: 



<9x 



a = CB 0fi 

= B 1)0 + B 0>1 , (A.10) 



/ Ft \ M / Ft x M+1 

= B N+ i^ M + Bn,M+U 

(A.ll) 



which allow us to obtain 



£ 2 6 = C 2 B 0)0 = CB lj0 + £-Bo,i 

= -^2,0 + + _B ,2, 

£ 3 5 = C 3 B 00 = CB 2 fi + 2CBii + CB Q2 

= £? 3 ,o + 3S 2 ,i + 3Si, a + ^,3, (A. 12) 

in 

" / ^m— n ! n J - > m— n,n- 

n=0 
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Here, the coefficients s % 
triangle, and thus 



(m = 0, 1, 2 . . . and n = 0,1, ... ,m) form a so-called Pascal's 



ml 



n\(m — n) 

From Eqs. ([A~9jl . ([A~T2jl . and (fA~T3|l . we have 



n=0 



A 2 "m!(£ m - n ■ 1) <9 n 6 
n! (m — n)\ dx n 



(A.13) 



(A.14) 



Substituting Eq. (IA.14D into Eq. (IA.6j) yields 



?Ti=0 



n=0 



n! (m — n)! 9x n 



(A.15) 



Using Eq. (|A.5|) . (C m n ■ 1) in this equation can easily be rewritten into (x m n ), which can 
further be rewritten as 



[X' 1 



(m-n)\ /<9 n x m 
ml V dx n 



(A.16) 



Substituting this into Eq. (|A.15|) yields 



m=0 



^ A 2n fd n x m \ d n b 



Moreover, because 



<9 n x m 
dx n 



n=0 



for n > m, 



dx n 



(A.17) 



(A.18) 



the summation over n = 1,2, ...,m in Eq. ()A.17|) can be extended to that over n 
1, 2, . . . , oo to obtain 



(ab) = J2 



m=0 



A 5 



n=0 



<9 n x m \ d n b 



dx n J dx n 



(A.19) 



Based on the commutativity between differentiations and filtering, after some mathematical 
operations we finally obtain Eq. (jA.l|) . 
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FIG. 1: Vector plot of the flow field described by Eqs. and (|7H|) in arbitrary units; note that 
this flow field is self-similar with respect to constant multiplication of the coordinates, (xi,X2) — > 
(0x1,0x2), where c is a real constant. The center point of the coordinate system shows the origin 
(x u x 2 ) = (0,0). 
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FIG. 2: Physical meaning of the negative diffusivity in the pure advection of a discontinuous step 
function in a stretching field. The solid lines in the upper figure denote fit = 0) and fit = t± > 0), 
and the dashed curves denote fit = 0) and f(t = t\ > 0). If / is advanced using a pure advection 
equation, its characteristic width gradually expands, as shown by the dots. A negative diffusion 
term must appear in the filtered advection equation to counteract this artificial expansion. 
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